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Summary 

The strain formulation in elasticity and the compatibility 
condition in structural mechanics have neither been under- 
stood nor have they been utilized. This shortcoming pre- 
vented the formulation of a direct method to calculate stress 
and strain, which are currently obtained indirectly by differ- 
entiating the displacement. We have researched and under- 
stood the compatibility condition for linear problems in 
elasticity and in finite element structural analysis. This has 
led to the completion of the “method of force” with stress (or 
stress resultant) as the primary unknown. The method in 
elasticity is referred to as the completed Beltrami-Michell 
formulation (CBMF), and it is the integrated force method 
(IFM) in the finite element analysis. The dual integrated 
force method (IFMD) with displacement as the primary 
unknown has been formulated. Both the IFM and IFMD 
produce identical responses. The IFMD can utilize the 
equation solver of the traditional stiffness method. The 
variational derivation of the CBMF yielded the existing sets 
of elasticity equations along with the new boundary com- 
patibility conditions, which were missed since the time of 
Saint- Venant, who formulated the field equations about 
1860. The CBMF, which can be used to solve stress, dis- 
placement, and mixed boundary value problems, has elimi- 
nated the restriction of the classical method that was 
applicable only to stress boundary value problem. The IFM 
in structures produced high-fidelity response even with a 
modest finite element model. Because structural design is 
stress driven, the IFM has influenced it considerably. A fully 
utilized design method for strength and stiffness limitation 
has been developed via the IFM analysis tool, and its merits 
and limitations are discussed. The method has identified the 
singularity condition in structural optimization and furnished 
a strategy that alleviated the limitation and reduced substan- 
tially the computation time to reach the optimum solution. 
The CBMF and IFM tensorial approaches are robust formu- 
lations because both methods simultaneously emphasize the 
equilibrium equation and the compatibility condition. The 
vectorial displacement method emphasized the equilibrium, 
while the compatibility condition became the basis of the 
scalar stress-function approach. The tensorial approach can 


be transformed to obtain the vector and the scalar methods, 
but the reverse course cannot be followed. The tensorial 
approach outperformed other methods as expected. This 
report introduces the new concepts in elasticity, in finite 
element analysis, and in design optimization with numerical 
illustrations. 


1.0 Introduction 

The theory of solid mechanics remained incomplete for 
more than one century. The deficiency pertained to the strain 
formulation in elasticity and the compatibility condition in 
structural mechanics. Our research has alleviated the defi- 
ciency. It is now complete for linear elastic problems. A 
veteran researcher should not be surprised over a deficiency in 
the strain formulation because some of the formulae and 
equations of the solid mechanics discipline were not com- 
pleted in the first attempt, but were perfected eventually. For 
example, perfecting the flexure formulae required more than 
one century between Galileo, Bernoulli, and Coulomb (ref. 1). 
Saint- Venant completed the shear stress formula that was 
initiated by Navier. Cauchy formulated the stress equilibrium 
equation that was also attempted by Navier in terms of dis- 
placement, but it contained only a single material constant 
instead of two (the rari-constant theory with one quarter for 
Poisson’s ratio). Green (ref. 2) subsequently resolved the 
misconception. We have completed the strain formulation 
(ref. 3) that was initiated by Saint- Venant in 1860. 

The compatibility condition (CC) is a significant ingredient 
in the theory of solid mechanics. Without a CC, the disci- 
pline would degenerate into a few applied mathematics 
courses in the analysis of determinate systems. The CC 
makes solid mechanics a research discipline that is practiced 
in the academia and in major research institutions throughout 
the world. The seed of the compatibility condition is 
ingrained in Hooke’s law, which forms the backbone of the 
theory of solid mechanics. The material law links stress {a} 
and strain {s} through the material matrix [k]: 

{G} = [K]{8} (1) 
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Stress must satisfy the equilibrium equation (EE). Like- 
wise, strain must satisfy the compatibility condition. That is, 
both stress and strain must simultaneously satisfy EE and 
CC, respectively, because one (stress) can be considered as a 
scaled version of the other (strain) and vice versa. Because 
stress is a tensor quantity (with six components), its calcula- 
tion requires an additional three compatibility conditions in 
addition to the three equilibrium equations. The tensorial 
nature does not allow the calculation of stress only from 
equilibrium considerations, and stress became an indetermi- 
nate quantity. The concept of indeterminacy is applicable to 
an elastic continuum as well as to a skeletal structure and a 
finite element model. A problem, either in continuum elastic- 
ity or in discrete structural analysis, is indeterminate if the 
number of stress unknowns n exceed the displacements m. 
For example, a plane stress elasticity problem is one-degree 
indeterminate because there are three stresses, n = 3, and 
two displacements, m = 2, with a one-degree redundancy, 
r = n - m = 1 . Stress determination can be represented by a 
symbolic equation that is obtained by coupling the equilib- 
rium equation and the compatibility condition: 



Figure 1. — Equilibrium equations and compatibility 
conditions in elasticity. 


Equilibrium Equation 
Compatibility Condition 


{Stress} 


J Mechanical Load 
[ Initial Deformation 


( 2 ) 


Equation (2) provides both necessary and sufficient condi- 
tion to calculate a high-fidelity stress state. Stress, calculated 
by any means or method, must be qualified through a com- 
pliance of equation (2); otherwise, such a stress state is liable 
to be incorrect. That is, the fidelity of stress determined via 
the stiffness method (ref. 4) has to be verified through a 
satisfaction of equation (2). Stress could not be calculated 
directly using equation (2) because the CC was not fully 
understood either in structures or in elasticity. 

We have studied the CC in structures (ref. 5) as well as in 
elasticity (ref. 6). The theory of solid mechanics is completed 
by adding the new information on the CC to the existing 
set of equations. The pie diagram in figure 1 graphically 
depicts the stress formulation and the strain formulation. 
Cauchy (1789-1857) developed the stress formulation 
during the first quarter of the 19th century. It contained both 
the field equations and the boundary conditions for an elastic 
continuum. Decades later, Saint-Venant (1796-1886) pro- 
vided the strain formulation (ref. 7) but only in the field of 
the continuum (shown as the third quarter in the pie chart). 
The boundary portion (marked as the shaded fourth quarter 
in the pie diagram), identified as the boundary compatibility 
condition (BCC), was missed until its formulation in 1986 
(ref. 3). The compatibility concept in structures, developed 
through a “cut” and “close” gap technique, is not even 


parallel to the available strain formulation in elasticity. The 
understanding and use of some aspects of the CC, which 
remained immature, is depicted as the shaded fourth quarter 
in the pie chart. 

Traditionally, solution to an indeterminate problem is 
obtained by utilizing the available information, contained 
within the unshaded three-quarter portion of the pie chart. 
Such approximate solution can be improved when informa- 
tion contained within the area of the entire pie diagram 
including the shaded area is utilized. The authors have 
formulated, verified, and utilized the additional information 
regarding the CC to solve problems. The available analysis 
methods for a solid mechanics problem are listed in table I. 
All analysis methods, including the flexibility method and 
the stiffness method, must satisfy both the EE and the CC; 
otherwise, the legitimacy of the predicted response can be 
questioned. 

We have improved the direct stress calculation method 
and the associated variational formulation in the theory of 
solid mechanics. In elasticity, the new method (ref. 8) is 
referred to as the “completed Beltrami-Michell formulation” 
(CBMF), with stress as the primary unknown (see top row in 
table I). The CBMF is applicable to both stress and dis- 
placement boundary value problems. It has alleviated the 
limitation of the classical Beltrami-Michell formulation that 
was only applicable to the stress boundary value problem. 
The method in finite element structural analysis is called the 
integrated force method (IFM) with force as its primary 
unknown (ref. 9). The IFM is an improvement over the 
classical flexibility method because the redundant concept of 
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TABLE I.— METHODS OF STRUCTURAL MECHANICS 
AND ASSOCIATED VARIATIONAL FUNCTIONALS 


Method 

Primary variables 

Variational 

functional 

Remarks 

Elasticity 

Structures 

Elasticity 

Structures 

Completed 

Beltrami-Michell 

formulation 

(CBMF) 

Integrated force 
method (IFM) 

Stresses 

Forces 

IFM variational 
functional 

Equilibrium and 
compatibility 
enforced 

Airy formulation 

Redundant 
force method 

Stress function 

Redundants 

Complementary 

energy 

Field 

compatibility 

enforced 

Navier formula- 
tion 

Stiffness 

method 

Displacements 

Deflections 

Potential energy 

Boundary 

compatibility 

noncompliant 

Hybrid method 

Reissner 

method 

Stresses and 
displacements 

Forces and 
deflections 

Reissner func- 
tional 

Boundary 

compatibility 

noncompliant 

Total formulation 

Washizu 

method 

Stresses, 
strains, and 
displacements 

Forces, defor- 
mations, 
and deflections 

Washizu func- 
tional 

Boundary 

compatibility 

noncompliant 


“cut” and “close” gap is not used, and it is applicable to 
dynamic analysis (ref. 10). Realizing the prominence of the 
stiffness method, a dual formulation referred to as the dual 
integrated force method (IFMD) (ref. 11) has been obtained. 
The equation structure of the IFMD and the stiffness method 
are similar in form and size, only the numerical values differ 
in the matrix [ D ] of the IFMD ([D]{X} = {/>}) and [K\ of 
stiffness method ([K]{X] = {/>}); here {X} and are 
displacement and load, respectively. The variational formu- 
lation of the new methods is referred to as the IFM varia- 
tional functional (n s ) (ref. 3). The stationary condition of the 
new functional yields both stress and strain formulations, 
new boundary compatibility conditions, and the displace- 
ment continuity conditions. 

The CBMF and IFM methods listed in the first row in 
table I can be considered as tensorial approaches in analysis 
because their primary unknown is either stress or internal 
force, which is a stress resultant. The tensor approach 
emphasizes both equilibrium and compatibility concepts. 
The stiffness and Navier’s displacement methods are vecto- 
rial approaches, which emphasize the equilibrium equation. 
The stress function approach is the scalar method, and the 
compatibility condition is its basis. The tensor approach can 
be specialized to obtain the vector and scalar methods, but 
the reverse course cannot be followed. For example, the 
equations of the new methods (IFM and CBMF) can be 
specialized to obtain the existing formulations of solid 
mechanics, such as Airy’s formulation, the stiffness method, 
the hybrid method, and the total formulation, as listed in 
table L Navier’s formulation cannot be transformed to obtain 
the CBMF in elasticity. Likewise, the IFM variational 


functional k s can be specialized to obtain the potential 
energy and complementary energy functionals, but the 
reverse course cannot be followed. 

Because structural design is stress driven, the IFM has 
become a useful tool in traditional fully stressed design as 
well as in design optimization. As it turns out, the IFM was 
devised to study the feasibility of full stress design in an 
indeterminate truss (ref. 12). An indeterminate truss cannot 
be fully stressed (ref. 13). This paper introduces the new 
compatibility concepts in elasticity, finite element analysis, 
and design optimization. This is followed by a conclusion 
section. A list of symbols and acronyms used in this report is 
given in appendix A to aid the reader, and appendix B 
contains a bibliography on the method of force for an 
in-depth study by an interested reader. 

2.0 Elasticity 

The material law is the foundation upon which the theory 
of solid mechanics is built. This stress (a)-strain (c) law 
(eq. (1)) was formulated and interpreted by Hooke (1635- 
1703), Young (1773-1829), Poisson (1781-1840), and 
others. The material law contained the essence of structural 
analysis. The constraint imposed on stress became the stress 
formulation, or the EE. Likewise, the condition on the strain 
became the strain formulation, or the CC. The stress and 
strain formulations along with the material matrix [k] given 
in equation (1) are sufficient for the determination of the 
stress state in an elastic continuum. The stress formulation is 
credited to Cauchy (1789-1857), who has formulated the 
field equation and the boundary condition, also known as the 
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traction condition. Saint-Venant (1797-1886) developed the 
strain formulation, and it was rederived by Sokolnikoff in its 
popular form (ref. 7). It was a general belief that the funda- 
mental elasticity relations were known for over a century. 
The thrust, therefore, was to develop approximate techniques 
because a closed-form solution can seldom be generated for 
the vast majority of the elasticity problems. In elasticity, 
Airy’s stress function technique (ref. 14) and Ritz’s dis- 
placement formulation (ref. 15) are two popular methods. 
Airy’s method is based on the compatibility concept, 
whereas Ritz’s method utilized the equilibrium equation. 
Their equivalents in the theory of structure and finite ele- 
ment analysis became the flexibility and the stiffness meth- 
ods, respectively. 

We now discuss the incompleteness in the strain formula- 
tion, which is quite straightforward to comprehend. Consider 
Cauchy’s stress formulation. It is complete because the 
treatment included the field as well as the boundary of an 
elastic continuum. In the field, 

x ij,j +bj=0 (3) 


BCC = 

( r)? r)? \ 

n r + 

(ds x 


[ dx dy y 

X 

v dy dx y 


where n x and n y are the direction cosines of the outward 
normal. With the BCC, the stress and strain formulations 
given above have parallel forms, each containing field equa- 
tions (eqs. (3) and (5)), and boundary conditions (eqs. (4) and 
(6)). The addition of the new BCC has completed the 
Beltrami-Michell formulation in elasticity. The completed 
method can be used to solve all three types of boundary value 
problems in elasticity: the stress boundary value problem, the 
displacement boundary value problem, and the mixed bound- 
ary value problem. The triviality of the known strain formula- 
tion in the field (f cc (u,v) = ,f{uy) - = 0, see eq. (5)) 

justified Navier’s method, which is obtained by rewriting 
Cauchy’s equilibrium equation in the displacement variables. 
Navier’s method can satisfy the compatibility condition in 
the field but not on the boundary because the new BCC is 
not automatically satisfied when expressed in displacements 
(; u and v): 


where Xy is the stress tensor and b t is body force. On the 
boundary, 


x ij n j = Pi ( 4 ) 

where rij is the direction cosine of the outward normal and p t 
is a traction component. 

Likewise, the strain formulation should have contained 
conditions in the field as well as on the boundary. However, 
it was formulated only in the field. The strain formulation in 
the field f cc for a plane elasticity problem can be written as 


BCC = 


a 2 v 

dxdy 


1 8 (du dv ' 

2 dyydy dx y 


d 2 u 

1 d 

' dv du^i 

dxdy 

2 dx ' 

H 

1 


( 7 ) 


The consequence of this noncompliance has to be investi- 
gated. The usefulness of the BCC is further investigated 
considering the CBMF in polar coordinates because this is a 
popular topic in the theory of elasticity. 


, d \ d ^i o 

dx 2 dy 2 dxdy 


( 5 ) 


where z x , s y , and y xy = 2z xy are strain components. Comparing 
stress and strain formulations we observe 

(1) Field equations are available for both stress and strain 
formulations (eqs. (3) and (5)). 

(2) The boundary condition is available only for the stress 
formulation (eq. (4)) but not for the strain formulation. 

The strain formulation remained incomplete since the time 
of Saint-Venant. That is, the theory of elasticity was incom- 
plete since 1860. Because of this deficiency, the classical 
Beltrami-Michell formulation could not be used to solve an 
elastic continuum with a displacement boundary condition. 
The strain formulation on the boundary has recently been 
completed (ref. 3). The boundary compatibility condition 
(BCC) for the plane stress problem can be written as 


2.1 CBMF in Polar Coordinates 

The strain formulation (that includes the field and bound- 
ary CCs) is required in the theory of elasticity because both 
stress and strain are tensor quantities. A variational approach 
has to be followed to derive the strain formulation, along 
with other equations for a two-dimensional continuum in 
polar coordinates because the boundary compatibility condi- 
tion cannot be generated from a nonvariational direct ap- 
proach. The stationary condition of the IFM functional yields 
both the stress and the strain formulations as well as the 
displacement continuity condition — or all equations of the 
CBMF. Consider a plane elastic continuum with uniform 
thickness t in the (r, 0) polar coordinate system, as shown in 
figure 2. It is made of an isotropic material with Young’s 
modulus E and Poisson’s ratio u. It is subjected to external 
load ( P r and P § ) along the boundary segment On the 
remainder of the boundary (f 2 )> the displacements are 
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Figure 2. — Two-dimensional elastic continuum 
in polar coordinates with loads P r and Pq 
and prescribed displacements u and v. 


restrained ( u = 0 and v = 0 ). Three stresses (c> r , a e , x r0 = x) 
and two displacements ( u and v) are the unknowns of the 
problem. Three strains (s r , s e , y r0 = y) can be calculated by 
scaling the stresses using Hooke’s law. The stationary 
condition of the functional n s of the IFM (ref. 6) with respect 
to displacement and stress function yields both the stress and 
strain formulations as well as the displacement continuity 
condition. The functional has three terms A, B , and W as 
follows: 


n s = a + B-W 


( 8 ) 


B = 


-§> 


^0s r d 2 (p ®s r <3(p ®sq5 2 (p 


D 


r 2 d0 2 r d r 


dr z 


^ y5 (p M y5 2 (p 


trdrdQ 




r 2 50 

rdrdQ J 



r 

fP _ 

m -) 

r 

(m_ ^ 

w = 

T 

liPy. + 

vP e 

d/:'! + (t 

) uR r + vRq 


_ A 


J 

^2 

v y 


^)(b r u + Z?0v)(rdrd0) 


D 


(9c) 


The domain ( D ) has boundary £ = £ i + 1 2 • Body forces 
are ( b r and Z? 0 ). Along the boundary segment £\, loads P r 
and are prescribed, whereas displacements ( u and v) are 
set free. The segment £ 2 has prescribed displacements u 

and v that can induce reactions R r and Rq. The derivation 
sets the uniform thickness to unity (t = 1) without any 
consequence. Body force potential V is defined as 

r dv dv 

b r = and b r = (10) 

dr rdQ 


The stress function cp is defined as 

rdr r 2 dQ 2 


(11a) 


The term A (a, u) represents the strain energy, expressed in 
terms of stress a and displacement u. The term i?(c, cp) is the 
complementary strain energy written in strain c and the 
stress function cp. The potential of the work done is W. The 
definitions of A, B , and W follow: 


d 2 cp 


-V 


T = ■ 


d 

dr 


f dcp ^ 

JdQ 


(lib) 


(He) 


A = 


0 a r du S x du 

h h 

dr rdQ 

0 x dv 0 erg dv 0 ti 


§1 

D 


0 


dr 


rdQ 


GqU 


trdrdQ 


J 


Details of the derivation are given in reference 6 and are 
not repeated here. The final expression of the stationary 
condition of the functional with respect to displacement and 
stress function has the following form: 
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8n s = (£j) [(field EE ) 8 [u } + (field CC) 8 (p]ds 
D 

+ (£) [(boundary EE ) 8 { u } ]df 

k 

+ ^[(boundary CC) 8 cp]df 

t 

+ (£) (continuity condition) 8 {reaction} dl = 0 


( 12 ) 


Field compatibility condition: 


1 d 2 a r 
r 2 ae 2 


d 2 a r 

dr 2 


(l + 2 o) do r d 2 


7 0 


dr d r 2 


o d a 0 ^(l + 2 o)^a 0 (l + u) d 2 x 


r 2 ae 2 
(l + u) dx 
2 dr 


dr 


drdQ 


= 0 


( 15 ) 


The field EE and field CC are the coefficients of the varia- 
tional displacement and stress function, respectively, in the 
surface integral terms of the functional in equation ( 12 ). 
Likewise the boundary EE (or traction condition) and the 
new boundary CC are the coefficients of the variational 
displacement and stress function, respectively, in the line 
integral terms. The continuity conditions are the coefficients 
of the variational reactions. The CC becomes relevant when 
the variational stress function is nonzero (Sep ^ 0) or the 
domain is indeterminate. The CC should not be enforced on 
a determinate domain with Sep = 0. 

The field equations and boundary conditions of the CBMF 
recovered from the stationary condition of the variational 
functional are listed below. The CCs recovered in equa- 
tion ( 12 ) are in strains, which are transformed into stress 
using Hooke’s law. 


Boundary compatibility condition: 


d / \ (l + u ) f 

-(ae-ua r )-— --a e+ a, 


/, / ,da e 2 ( 1 + u) 

v ’ dr K ’ 50 r 


«e = 0 


The kinematics boundary conditions are 


u-u = 0 
v-v =0 


where 


(16) 


(17a) 

(17b) 


Completed Beltrami-Michell Formulation in Polar 
Coordinates: 


Field equilibrium equations: 


dr r <90 r 


dx 1 doQ 
dr r <9© 


+ — + Z? 0 = 0 


Boundary traction conditions: 

n r a r +n§ x = P r 


(13a) 

(13b) 


(14a) 


r and 0 

polar coordinates 

a r , ce, and x 

stress components 

s r , £e, and y 

strain components 

u and v 

displacements 

b r and be 

body forces 

P r and Pq 

tractions applied along boundary segment 1 1 

u and v 

initial displacements along boundary segment 

n r and n$ 

z 2 

direction cosines of the outward normal 


2.1.1 Verification of boundary compatibility condition . — 
Green’s theorem is used for a quick verification of the new 
BCC. The BCC is inserted in the line integral coefficient to 
recover the well-known field CC in the surface integral term. 
The integral theorem in polar coordinates can be written as 


n r x + « e a e =P e 


(14b) 



d g e 

rdQ 


ds = j)(G r « r + GqHq )d l 


(18a) 
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where G r and Gq are the coefficients of direction cosines n r 
and riQ, respectively. The coefficients are defined for the 
CBMF as 




Figure 3. — Annular plate subjected to 
linear temperature distribution. 

2.1.2 Illustrative example: annular plate subjected to 
thermal load . — We will now illustrate the CBMF calcula- 
tion strategy and the use of the BCC through the solution of 
a radially symmetrical annular plate with mixed boundary 
conditions for thermal load. Consider a plate made of an 
isotropic material with Young’s modulus E and Poisson’s 
ratio u. It has thickness t (considered unity) with outer and 
inner radii of a and b , respectively, as shown in figure 3. The 
inner boundary is fully restrained while the outer boundary is 
free to expand. The annular plate is subjected to a tempera- 
ture distribution: 


(18f) 

The right-hand coefficient in equation (18f) is the field 
CC. The compatibility concept applies to the field as well as 
to the boundary. The analytical form of the compatibility 
expression changes in compliance with the domain and the 
boundary. The same interpretation is true for Cauchy’s field 
EEs. 

The BCC written in displacements is not a trivial 
condition: 


T = T b + ^ Ta ^ (r-b) (20) 

\ a ~ b ) 

The temperature has a linear variation with values T a and 
T b at r = a and r = b , respectively. The coefficient of thermal 
expansion is a. The CBMF equations for the annular plate 
subjected to a thermal load are given below: 

Equilibrium equations 


d 2 u 

r d2y 

dv 

59 2 

r drd% 

50 


d 2 u 2 d 2 v dv du 

+ r h v-r — — - r r — 

dr 50 dr 50 

V ur J 


n Q =0 


(19) 


Field: 

aa r | (a r -cje) q 
dr r 

Traction or boundary: 


(21a) 


The field CC is a second-order differential equation in 
strains, while the boundary counterpart is a first-order 
equation. This characteristic is applicable to the stress 
formulation. The field EEs are first-order equations in stress, 
while the boundary (or traction) conditions are algebraic 
equations. The BCC is not a continuity condition in dis- 
placement. The BCC is an independent condition. It forms a 
new elasticity expression that was missed since the time of 
Saint- Venant. 


a r = 0 at r = a (21b) 

Compatibility conditions 

Field: 

d , v (l + o j , v d T 

— (CTe-UCT r ) + (CT e -q r ) = -a£'— (22a) 

dr r dr 
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Boundary: 


G0 - oa r = -a ET at r -b (22b) 

Displacement boundary condition 

u = 0 at r = b (23) 

The temperature term is contained explicitly in both field 
and boundary CCs (see eqs. (22a) and (22b)). The EE and 
CC in the field are rearranged to obtain the following two 
simple working equations: 

d , v dr ^ x 

(a r +ae) = -aE— (24a) 

dr dr 


da r | (^r-^e) Q 
dr r 


(24b) 


The CBMF solution strategy has two distinct steps. First, 
the stress state is calculated from the field EE and CC given 
by equations (24a) and (24b), for the boundary conditions 
given by equations (21b) and (22b). The displacement 
boundary condition (given by eq. (23)) is not used to calcu- 
late the stress state. Pure displacement cannot induce stress. 
Only the derivative of displacement, which is strain, pro- 
duces stress. The stress solution follows: 


fa(fl-r) 

3 r 2 (a-b) a 2 + b 2 +o^ 2 -b 2 j 

r 2 L 2 + b 2 + u(a 2 -b 2 )]' 

x T a < +rb 2 [a (l-o)-Z?(2-o)] > (25a) 

+ab 2 [a(l-o)-Z?(2-o)J 

V l 

+T b < +rb (2 + o) — b{\ + oj 

-\~ub b (l + o) + u (2 + 



E a 

3r 2 [b-a ) a 2 +b 2 +u[a 2 -b 2 j 


f 

2r 3 

a 2 +b 2 + u^a 2 -b 2 j 

Ta' 

-r 2 

(? +2 b 3 + o|a 3 -6 3 J 

V 

+a 2 i 

b 2 [a(l-u)-b(2-u)] 


-2 r 3 a 2 +b 2 +u(^a 2 -b 2 j 
+ T b < +r 2 a 3 -b 3 +3ab 2 + o|<2 3 -b 3 j 
+a 2 b 2 [a (2 + o)-Z?(l + o)] 


(25b) 


Next, the displacement function is back-calculated by 
integrating the stress. The displacement boundary condition 
u = 0 at r = b is used to evaluate the integration constant. The 
displacement function is 

a(l + u)(Z?-r) 

3 r 2 (b-a) a 2 +b 2 +v(^a 2 -b 2 j 

r 2 a 2 +b 2 +o|a 2 -b 2 j 

x T a < +ra 2 |^a-2Z?-o(a-Z?)J ► (26) 

+a 2 Z?[a(l-o)-Z?(2-o)] 

V l 

-r 2 a 2 +b 2 + 0 ^a 2 -b 2 j 

+ T b < +ra 2 [2 a-6 + o(a-fe)J 
+a 2 Z?[a(2 + o)-Z?(l-Mj)] 



The numerical values of the response parameters for T a = 
100 °C, T b = 50 °C, and a = 12xlO“ 6 /°C are 

(1) at r = a: o r = 0 ksi, g 9 = -17.5 ksi, and u = 0.012 in. 

(2) at r=b: <j r = 14.2 ksi, a e = -13.7 ksi, and u = 0 in. 
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The sum of the stresses a r + a e = 18508 - 1800r has a 
linear variation with respect to the r-coordinate because of a 
similar distribution of temperature (see eq. (20)). The solu- 
tion steps in the CBMF are discussed in appendix C. 

2.2 Treatment of Boundary Compatibility in Na- 
vier’s Displacement Method 

The displacement method of Navier is universally ac- 
cepted, and it has been applied to solve many problems. It 
has field equations and an appropriate set of boundary 
conditions. Adding an extra BCC (see eq. (19) for example) 
may appear to hamper the method with an excessive number 
of boundary conditions for which no analytical solution can 
be generated. This observation is an aberration but has no 
serious implication. We discuss the BCC by solving a plate 
flexure problem because this is a popular example and its 
solution occupies a good portion of textbooks (ref. 16). The 
rectangular plate is made of an isotropic material with 
Young’s modulus E and Poisson’s ratio o, thickness t , and 
spans of 2 a and 2c along x- and y-coordinate directions, 
respectively, as shown in figure 4. It is subjected to a distrib- 
uted transverse load q with a sinusoidal variation along the x- 
coordinate direction, while being uniform along the y- 
direction as 


q = qo cos(nx/2a ) 


(27) 


EE: 


a 2 M v d 2 M v d 2 M : 


xy 


dx z 


dy 1 dxdy 


(28a) 


CCl : — (M v -uM x )-(l + u) — M xy = 0 (28b) 


dx 


dy 


CC 2 : — (M x - uM v ) - (1 + u) —M xy =0 (28c) 


dy 


dx 


The three field equations are reduced to two by eliminat- 
ing M xy between the two CCs. The working equations of the 
CBMF and the boundary conditions are listed below: 

Field equations of the plate flexure problem in CBMF 


V 2 (M x +M y ) = {\ + »)q 


(29a) 


S 2 \ d 2 


—(M y - u M x ) ~(m x - vM y ) = 0 (29b) 

dx dy 


Boundary condition of the plate flexure problem in 
CBMF 


The plate is simply supported in the x-direction along AD 
and BC, while it is clamped in the y-direction along AB and 
DC. The three moments M x , M y , and and the single 
transverse displacement w are its response variables. The 
CBMF formulation for the plate flexure problem is obtained 
from the stationary condition of the IFM functional (ref. 3). 
The stationary condition yields three field equations in 
moment variables, containing one EE and two CCs: 



Figure 4. — Rectangular plate in flexure. 


Simply supported edges (AD and BC): 


M x = 0 (30a) 

M y = 0 (30b) 

Clamped edges (AB and DC): 

M x - vMy = 0 (31a) 

M xy = 0 (31b) 


The field equation and boundary conditions of Navier’ s 
method, expressed in displacement (w), are as follows 

Field equation in Navier’ s (or stiffness) method 

V 4 w = — (32) 

D 

No explicit field CC because it is automatically satisfied 
in displacement. 
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Boundary condition in Navier’s (or stiffness) method 

Simply supported edges (AD and BC): 


d 2 w d 2 w 


dx z 


- + u- 


ay 2 


= 0 


w = 0 

Clamped edges (AB and DC): 

^ = 0 

ay 

w = 0 


(33) 


(34) 


Boundary condition in CBMF written in displacement 

Simply supported edges (AD and BC): 


a 2 w 

ax 2 

a 2 w 

ay 2 




= 0 


(35) 


Clamped edges (AB and DC): 


a 2 w 

ax 2 

d 2 w 

dxdy 


= 0 


= 0 


(36) 


In the CBMF, the kinematics condition w = 0 is not used 
to calculate the moments because it represents a rigid motion 
and does not induce moment or stress in the plate. During 
the back-calculation of displacement from moments, the 
CBMF method uses the rigid body displacement (w = 0) 
along the boundaries AB, BC, CD, and DA. The slope 
condition is not at all used in the CBMF. The field equations 
in the CBMF and the Navier’s method are equivalent to a 
fourth-order equation. Each method has two boundary 
conditions at each edge. The boundary conditions are consis- 
tent with the order of the field equations for both methods. In 
other words, the new boundary compatibility conditions do 
not create an infeasibility issue. 

2.2.1 Solution strategy . — In the CBMF the second-order 
differential equation (29a) is solved first to obtain the sum of 
moments M x + M y . Then equation (29b) is uncoupled and 


solved independently. The CBMF requires the solution of 
two second-order uncoupled differential equations. It is 
simpler than the stiffness method, where a fourth-order 
differential equation must be solved. The CBMF solution is 
obtained using separation of variables. 


M x -M x (y) cos ax 

(37a) 

My -My (y) cos ax 

(37b) 

M xy = M xy (j) sin ax 

(37c) 


The CBMF moment (M x , M y , and M xy ) solution follows: 


M x - cos ax 


C3 {e ay + e~ ay j 


-a 


( 1 - 0 ) 

2(1 + u) 


c iy (e a y-e- a y)-^ 

' ' a 


M y = cos ax 


(q-c 3 )( 

_a-u 

2(1 + u) 


e ay +e -ay 


+a — — — Qj ^e a:y 


a 


, , sin ax 

M ,v= — 


(2C 3 —C\)^e ay -e~ ay ^ 




(38a) 


(38b) 

(38c) 


Displacement w is calculated by integrating a moment 
component. The integration constants are evaluated from the 
rigid-body boundary condition w = 0. The displacement 
function follows: 


- cos ax 


a 2 D 


H-) 


C 3 ( l + v)(e a y+e- a y 
-uC'i (c a >’ + e - a >'j 


(39a) 


-a 


e ay +e -ay 
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The constants C i and C 3 are given by 


Q= ^ "7 (1 + 0) 


e ac - e “ ac 
' e 2ac + 4 ac - e _2ac 


0 2 W 0 2 W 


(l + u) 


e ac — e~ ac j + a (l — u)c ye ac + e~ ac 
e 2ac + 4ac-e~ 2ac 


d 2 w _ d 2 w 
dx 2 dy 2 


’M x -M y = 0 


Identical solutions are obtained for displacement and 
moments by Navier’s method. 

2.2.2 Boundary conditions . — The moment boundary 
conditions in the CBMF are expressed in terms of curvatures 
using the moment curvature relations to obtain equations 

(35) and (36). The CBMF boundary conditions, when 
expressed in displacement, do not identically map into the 
conditions of the displacement method (see eqs. (33) to 

(36) ). In the CBMF, the boundary conditions are expressed 
in moments or curvatures, but in the displacement method 
these are written in moment, displacement, and slope. 
Consider for example, the simply supported boundary 
conditions along AD and BC. The CBMF method imposes 
zero-moment conditions M x = M y = 0 along the boundary. 
When expressed in displacement these become zero condi- 


tions m curvatures 


d 2 w d 2 w 


The boundary 


conditions in the CBMF and displacement methods appear to 
be different. The same logic applies to the clamped bound- 
ary. In the CBMF, conditions are expressed in curvatures 


= 0 , while they are in terms of displace- 


ment and slope in Navier’s method w = = 0 . How- 

{ dy ) 

ever, the response (M x , M y , M xy . , w) is the same by either 
formulation because the boundary conditions of the two 
methods can be shown to be equivalent. This is explained 
next. 

For simply supported boundaries along AD and BC (x = 
-a and x = a, and -c < y < c), the displacement can be 
differentiated with respect to the y-coordinate without any 
consequence because the displacement is uniform across this 
boundary. 


dw d 2 w 


w = 0^> — = 


dy dy 2 


For clamped boundaries along AB and DC (y = -c and y = 
c , and -a < x < a), we can differentiate the slope and dis- 
placement with respect to the x-coordinate. 

dw d 2 w 

— => = 0 (41a) 

dy dxdy 


2 

(iv = 0) => — =^-^- = 0 

dx dx 2 


Mjcy =0 


~VM y =0 


The boundary conditions become equivalent in Navier’s 
and the CBMF methods. Thus both methods produce an 
identical response. The equivalence concept may become 
invalid if the rigid (or uniformity) condition (used in 
eqs. (40) and (41)) is not allowed, such as along a flexible 
finite element interface (ref. 17). All boundary conditions, 
including the continuity conditions, should be examined for 
equivalence between the CBMF and displacement method. A 
comprehensive comparison is not included in this paper 
because it requires substantial discussion. 

2.3 Summary on Elasticity 

In the theory of linear elasticity, Navier’s method solved 
all three types of boundary value problems in elasticity, 
whereas the classical Beltrami-Michell formulation could not 
solve the second and third boundary value problems. In other 
words, the classical theory of elasticity was incomplete. 
Love’s statement below was applicable only to the stress 
boundary value problem (ref. 14): 
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“It is possible by taking account of these [compati- 
bility] relations to obtain a complete system of equa- 
tions [Beltrami-Michell Formulation] which must be 
satisfied by stress components, and thus the way is 
open for a direct determination of stress without the 
intermediate steps of forming and solving differential 
equations to determine the components of displace- 
ments.” 

A telltale sign pointing to the incompleteness was ignored. 
The BCC has alleviated the limitation of the classical 
Beltrami-Michell formulation. The completed method, 
CBMF, can solve all three classes of boundary value prob- 
lems in elasticity. However, we might have exposed a new 
issue: for example, for a plate flexure problem, the boundary 
condition should be imposed on moments in the CBMF, 
which translates to curvatures in displacement method 
instead of the popular slope condition. Resolution of the 
issue requires further research. 

3.0 Structures 

The traditional compatibility condition in structures was 
formulated through a concept of redundant force. The CC 
was improvised by cutting a redundant member and subse- 
quently closing the gap. The process automated by Denke 
(ref. 18) was referred to as “structure cutter.” The fictitious 
“cut” and “close” gap technique is not even parallel to the 
then available strain formulation in elasticity. The force 
method based on the concept of redundant forces has disap- 
peared for all practical purpose because it was clumsy and 
could not be extended to dynamic and buckling analyses 
(ref. 19). We have extended the strain formulation to struc- 
tures (ref. 20). The CCs thus obtained are added to the 
equilibrium equation. This method with force as the primary 
variable is referred to as the “integrated force method” (IFM) 
(ref. 5). The IFM bestows simultaneous emphasis on the EE 
and the CC. Conceptually, the CBMF and IFM are parallel 
formulations in elasticity and in structures, respectively. To 
utilize the vast resources that have been developed for the 
regular stiffness method, like equation solver as well as pre- 
and postprocessors, a dual method to the primal IFM has 
been formulated. The dual integrated force method (IFMD) 
(ref. 11) with displacement as the primary unknown is 
obtained by transforming the equations of the IFM. Govern- 
ing equations of the IFMD are symmetrical and resemble the 
stiffness equation in size and form. Both the IFM and IFMD 
produce an identical response. The IFMD can utilize the 
equation solvers of the popular stiffness method. In fact, a 
stiffness method code can be converted into IFMD software 
with minimal programming effort. 

The formulation of the CC and IFM are introduced here 
through a simple example. Then the equations of the IFM 
and IFMD are given. Fidelity of the IFM and IFMD solu- 
tions is discussed. 


3.1 Illustrative Example 

The basic steps of the IFM that included the CC are intro- 
duced, considering the example of the three-bar truss shown 
in figure 5. The truss is made of steel with Young’s modulus 
E = 30 000 ksi, bar areas A x = A 2 = 1.0 and A 3 = 2.0 in. 2 , and 
load P x = 50 and P y = 100 kip. The n = 3 bar forces F u F 2 , 
and F 3 and the m = 2 displacements X x and X 2 are its five 
response variables. The structure has m = 2 equilibrium 
equations, which are written at the free node 1 : 


1 

VT 

-i 

VT 


-i 


-i 

VT 

-i 

VT 


f 2 

f 3 



(42a) 


[5]{F} = {P} 


(42b) 


Here, [B] is the rectangular equilibrium matrix of dimen- 
sion mxn, {F} is the ^-component force vector, and {P} is 
the m-component load vector. Three forces cannot be deter- 
mined from the two EE. One CC is required. This is obtained 
as an extension of the strain formulation in elasticity to 
structures (ref. 20). The CC is generated in three steps. First, 
the deformation-displacement {P}-{A} relation (DDR) is 
obtained as (P) = [B] T {X}. Then displacements are elimi- 
nated from the DDR to obtain the CC in terms of deforma- 
tions. Finally, the deformations are eliminated in favor of bar 
forces using the flexibility relationship to obtain the CC in 
forces. 


100 4 * 100 

0 m a 



Figure 5. — Three-bar truss with loads P, bar forces F, and 
displacements Xi and X 2 . Dimensions are in inches. 
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The DDRs {P} = [B]T {X} that are analogous to the strain 
displacement relation in elasticity follow: 


V 

< p 2 > 

p 3 . 


1 

77 

o 

i_ 

"77 


'77 

-i 

i 

'77 


X\ 

X 2 


(42c) 


77 0 

1 -1 


-1 

VT 

-1 

a 

2 


A 


50 

Fl 

> = < 

100 

TV 


0 


(42i) 


Solution to equation (42i) yields the forces as 


Here, p b p 2 , and p 3 are the bar deformations, andA^ and 
X 2 are nodal displacements. The CC is obtained by eliminat- 
ing the two displacements from the three strain displacement 
relations: 


>r 


' 5.025 ' 

<f 2 

> = < 

42.893 > 

TV 


75.736^ 


[l -77 l] 


Pi 

PlHo} 


[P 3 

[C]{p} = {Stf} 


The displacements are back calculated from forces using a 
(42d) formula {X} = [J\[G] {T 7 }, where \J\ = m rows of [S]~ T . 


(42e) 


\u\ -fiopy f 2 i _ r°- i4 3i 

[vjl E ){2F l -F 2 J io.llOj in 


(42j) 


Here, [C] is the rxn (where r = n - m) compatibility 
matrix, and {87?} is the r-component initial deformation 
vector for the problem; for a mechanical load, {87?} =0. The 
CC ([C]{P} = 0) is analogous to the strain formulation in 
elasticity. 

The deformations are related to the forces {P} = [G]}/ 7 } 
through the flexibility matrix [ G ] of dimension nxn : 


pr 

" P2 ’ 

JA 



0 

0 


0 



-] 

AEj 3 


f i 

\<F 2 > 

TV 


(42f) 


The generalization of the calculations shown for the three- 
bar truss became the IFM. The IFM, which is the discrete 
analogue of the CBMF, simultaneously emphasized the EE 
and the CC. The IFM was not developed earlier because the 
CC was not available. 

3.2 Equations for the Integrated Force Method 

In the IFM of analysis, a finite element model, is charac- 
terized by two attributes: n forces {T 7 } and m displacements 
{X}. The IFM has one set of equations to calculate the 
internal forces and another set to back-calculate the dis- 
placements. The IFM equations for static and free-vibration 
(or frequency) analyses are as follows: 

Static analysis: 


The CC is expressed in forces by eliminating deformation 
between equations (42d) and (42f): 


[S]{F} = {/>*} (43) 


f 






(42g) 


or 

[C][G]{F}={5tf} = {0} (42h) 

The EE = \P\) and the CC ([C][G]{F} = {5 R) = 

{0}) are added to obtain the governing equation of the IFM 
as[S]{F} = {P}: 


W=F]{[ g ]T}+{p 0 }} ( 44 ) 


where [S] = 


[*] 


[j] = m rows of f[S] 1 j , and 


( * i 

(A 1 

7H 

II 

1 

o 
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Free- vibration analysis: 


Frequency analysis: 


[S]-co 2 


M'M 

“Tor - 


TM°} 


J J 


(45) 


Calculation of displacement mode shape in frequency 
analysis can use the displacement calculation formula given 
by equation (44) in static analysis. The symbols in the static 
and free-vibration analyses are listed as follows: 

{F} n component internal force vector 

{F} m component external load vector 

{8R} (r = n — m) component effective initial deformation 
vector 

{P°} n component initial deformation vector 

(P) n component deformation vector 

{X} m component displacement vector 

[Si nxn IFM governing matrix 

IB] mxn equilibrium matrix 

\Q rxn compatibility matrix 

[G] nxn flexibility matrix 

[J] mxn displacement coefficient matrix 
\M] mxm mass matrix 

m number of displacement unknowns 

n number of force unknowns 

co circular frequency 

3.2.1 Dual Integrated Force Method . — Solution of a 
symmetrical set of equations is very well developed, espe- 
cially because of its use in the finite element stiffness 
method. Furthermore, the Patran pre- and postprocessing 
software (MSC. Software Corporation) is a very useful 
analysis tool. The IFM can benefit from the solver and the 
software, provided the IFM equation can be cast in the 
format of the stiffness method. This has been achieved, and 
the method is referred to as the dual integrated force method 
(IFMD). The equations of the IFMD are obtained by map- 
ping forces into displacement (ref. 11). Like the IFM, the 
dual method also has two sets of equations. The first set is 
used to calculate displacement, while the second set back- 
calculates forces. The primal and dual methods produce 
identical solutions for force, displacement, and frequency. 
The equations of the IFMD follow: 


[Z)]-co 2 [M] {x} = {0} 


(47) 


Here, [D] is an mxm governing matrix of the IFMD, and 
{F 0 } is the initial load vector. Calculation of force mode 
shape in frequency analysis can use the force calculation 
formula in static IFM analysis given by equation (46b). The 
governing equation of the IFMD resembles that of the 
regular stiffness method, [K\{X] = {F stlff }, where [K] is the 
mxm stiffness matrix, and {F stlff } is the m-component load 
vector of the stiffness method. Both [K\ and [ D ] are symmet- 
rical matrices with identical dimensions, but their coeffi- 
cients are not equal (Dy ^ Ky). Notice the flexibility matrix 
[G], which is maintained in a pristine state in both the primal 
and the dual methods. This makes the calculation of design 
sensitivity straightforward via the IFM and IFMD. 


3.2.2 Other methods . — The matrices of the IFM can be 
used to formulate the mixed method and the total formula- 
tion of structural mechanics given in table I. The reverse 
course cannot be followed. For example, the equations of the 
stiffness method cannot be transformed to obtain the IFM. 
However, the IFM equations were specialized to obtain the 
IFMD, which in essence is the displacement method. Like- 
wise, the mixed method and the total formulation can be 
constructed from the IFM equations. 

3.2.2. 1 Mixed method : The mixed method treats force 
{F} and displacement { X) as the simultaneous unknowns. 
The governing equation of this method is obtained by 
concatenating the two IFM equations as 


[S] [0]1 


H ] 

L[-M[g]] [/]J 

wr 



(48) 


where [7] is the identity matrix. 

3. 2. 2. 2 Total formulation’. The total formulation consid- 
ers force {F}, displacement {X}, and deformation {P} as the 
simultaneous unknowns. The governing equation of the total 
formulation is obtained by adding the flexibility formula 
[G] {F} = {P} to the equations of the mixed method. 


Static analysis: 


[Z)]{x} = {/>} + {/> 0 } (46a) 

{^} = [G]- 1 {[5f{x}-{p°}} (46b) 

where [D] = [5][G] -1 [*f and {p 0 } = [5][G] _1 {p 0 } 


.Co 

o 

o 

T) 


( * ) 
F 

-F][G] [/] 0 

<{x} 

> = < 

M{e 0 } 

-[<?] 0 [/]] 

[{P}J 


{»! 


(49) 


The discussion to follow emphasizes only the force 
method (IFM) and the dual method (IFMD). The popular 
stiffness method is used to compare results. 
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3.2.3 Solution fidelity. — The stiffness method is a popular 
research topic throughout the world. Commercial code 
developers are continuously improving the performance of 
this method through techniques like reduced integration, the 
use of a bubble function, and others. Also, the IFMD concept 
might have been incorporated (ref. 21) into the stiffness 
method. It is difficult to numerically compare with such 
improved techniques in the commercial software because of 
too many attributes, some of which may even be classified. 
For a meaningful comparison in a controlled environment we 
have developed the IFM/ Analyzers software (ref. 22). This 
finite element code incorporates the IFM and IFMD as well 
as the standard stiffness method. The code has a total of 44 
different types of elements for each method. The generation 
of each element involves similar energy expressions as well 
as the same interpolation polynomials and numerical integra- 
tion technique. The IFM/ Analyzers element library includes 
beam, membrane, plate, and solid elements of different 
shapes. Some elements have midside nodes. IFM/ Analyzers 
is a Fortran 77 code, written for both sequential and parallel 
calculations. It performs linear analysis for mechanical and 
thermal loads and initial displacement as well as frequency 
analysis. The analyzer uses the NASA GPS sparse solver 
(ref. 23) and Harwell routines (ref. 24). Solutions have been 
obtained for finite element models up to one-half million 
unknowns. The IFM/ Analyzers code has also been reduced 
to obtain a modest code with all three methods. The small 
code named “IFM-UE” can be used for undergraduate 
education in engineering. The code accompanies the text- 
book, reference 5. The solution capacity of this code is about 
5000 equations. It has five different types of elements that 
can model skeletal frames and membrane structures. The 
IFM, IFMD, and stiffness method have been compared in a 
controlled environment for a set of test problems (refs. 11 
and 22). Solution is given for a few typical examples to 
illustrate the robustness of the IFM/ Analyzers code. 

3.2.3. 1 A cantilever plate with settling of support : A can- 
tilever plate ABCD is supported along AD with prescribed 
displacements u and v along this edge 


( 


*- 0 , --<_>><- 
2 2 


as shown in figure 


6(a). The pre- 


scribed displacement functions are 


u = 


4(2 + u) 
3Eh 2 


r 9 

v 4 


(50) 


v = 


4o L 9 
Eh 2 


(51) 


Mechanical load X applied along the free edge BC is 




Number of elements 


Figure 6. — Stress and displacement solutions for cantilever 
plate, (a) Cantilever plate, (b) Convergence of stress at 
location Q. (c) Convergence of displacement at location B. 

The example has been used to provide a benchmark solu- 
tion for finite element analysis. The analytical response, 
consisting of displacements u and v is as follows (ref. 25): 


T = 1-42— 


(52) 


4 

u = 1 

y{x 2 -2Lx)+ ^ + U y 

( h 2 \\ 

n 9 

— — y 

Eh 2 

3 

v 4 )_ 


(53) 
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Lx z 


Eh 1 


x 2 ( 

u v (x 

3 


4 + 5o 2 

■ L) H h x 

12 


(54) 


A finite element method solution was obtained using the 
IFM and stiffness method. A simple quadrilateral membrane 
element was used. The stiffness element had 8 displacement 
degrees of freedom. The IFM considered five force 
unknowns (ref. 5) with the displacement function that is 
identical to the stiffness method. The parameters in the 
problem were length L = 20 in., height h = 1.0 in., thickness t 
= 0.125 in., Young’s modulus E = 30 000 ksi, and Poisson’s 

ratio o = 0.3. Normal stress is calculated at the top fiber 

location Q with coordinates x = 1.0 in. and y = 0.5 in. as 
shown in figure 6(a). Transverse displacement v B is deter- 
mined at location B with coordinates x = 20 in. and y = - 
0.5 in. The analytical solution for stress and displacement are 
as follows: 


ct£ = 19psi 
v B =45.66 xlO -6 in. 

The finite element solutions are given in table II. Ten 
finite elements were used along the y-coordinate direction to 
accommodate the variation of load and prescribed displace- 
ment given in equations (50) to (52). Along the 
x-coordinate direction, the number of elements was progres- 
sively changed from 2 to 20. A 70-element model for exam- 
ple, had 7 and 10 elements along the length and depth, 
respectively. 


TABLE II. — STRESS AND DISPLACEMENT SOLUTIONS 
FOR CANTILEVER PLATE SHOWN IN FIGURE 6(a) 



Stress at Q, a,b 

psi 

Displacement v B at B, a c 
lO^in. 

Model 

IFM/IFMD 

Stiffness 

IFM/IFMD 

Stiffness 

2x10 

18.99 (100) 

5.18(27.0) 

43 (94.0) 

14 (30.2) 

4x10 

19(100) 

11.80 (62.0) 

45 (98.5) 

29 (62.9) 

6x10 

19(100) 

15.04 (79.0) 

45 (99.4) 

36 (79.1) 

8x10 

19(100) 

16.60 (87.0) 

46 (99.7) 

40 (87.0) 

10x10 

19(100) 

17.97 (94.6) 

46 (99.8) 

42 (91.2) 

20x10 

19(100) 

18.71 (98.5) 

46(100.0) 

45 (97.6) 


a Numbers in parenthesis represent percentages of analytical value. 
b Position on plate in figure 6(a), x = 1 andy = -0.5 in. 

Position on plate in figure 6(a), x = 20 andy = -0.5 in. 


The primal (IFM) and dual (IFMD) methods yielded iden- 
tical solutions. For the 20-element (2x10) model, the 

IFM/IFMD code yielded a stress value of =19 psi, 
which was identical to the analytical solution. The stiffness 
method, for the 2x10 model, yielded a stress of = 

5.2 psi, which was only 27 percent of the analytical solution. 
For other, dense models the IFM/IFMD result continued to 


coincide with the analytical solution as shown in figure 6(b), 
while the stiffness method converged to a stress value that 
retained a 1.5 percent residual error, even for a large model 
with 200 elements. For the 20-element model, the 
IFM/IFMD code yielded a displacement value of 43xl0“ 6 in. 
compared with the analytical solution of 46x1 0“ 6 in. The 
stiffness method, for that model, yielded a displacement 
value of v B = 14xl0“ 6 in., which was only 30 percent of the 
analytical solution. For other, dense models the IFM/IFMD 
result coincided with the analytical solution as shown in 
figure 6(c), while the stiffness method converged to a dis- 
placement value that retained a 6 percent residual error using 
the biggest model with 200 elements. 

3.23.2 A plate flexure problem : The rectangular plate 
under a concentrated load P = 1 kip, shown in figure 7(a), 
was a finite element test problem. It was made of steel with 
Young’s modulus E = 30 000 ksi, Poisson’s ratio u = 0.3, 
thickness t = 0.25 in., length a = 24 in., and width b = 12 in. 
It was clamped along all four edges. The analytical solution 
for the transverse displacement under the load was w P = 
2.42x1 0” 3 in., and the moment at x = 0, y = bl 2 was 
M b = 167 lbf-in./in. A four-node quadrilateral element 
(ref. 9) was used to solve the problem. Each node has 
3 degrees of freedom consisting of a displacement and two 
rotations (w h Q xi , 0 >7 ). A standard cubic polynomial with 12 
unknowns was used to approximate the displacement func- 
tion w(x,y). For the IFM/IFMD code, the force interpolation 
contained nine unknowns. Linear variation was considered 
for the normal moments M x and M y with four unknowns 
each, while the twisting moment was uniform across the 
element. The displacement and moment solutions for the 
plate are graphed in figures 7(b) and (c), respectively. 
Convergence was achieved for IFM/IFMD in displacement 
and moment for a model with eight elements. The stiffness 
method from two commercial codes were used. These 
methods required 64 or more elements for convergence with 
some residual errors. 

3.23.3 Frequency analysis of a turboprop blade : The 
turboprop blade, shown in figure 8(a), was made of steel 
with a weight density of p = 0.289 lbf/in. 3 It had a spatial 
distribution of mass with a uniform thickness of 0.1 in. It 
was about 10 in. long and 4 in. wide at the midspan location. 
An eight-node brick element was used to calculate its fre- 
quency. The element had 24 degrees of freedom that corre- 
spond to three displacements at a node. The IFM element 
(ref. 22) used a total of 18 force unknowns: 4 for each 
normal stress and 2 for each shear stress. The brick element 
of a commercial code was comparable to the stiffness 
method of the IFM/ Analyzers software. The frequency 
analysis results obtained by IFM/IFMD are given in table III. 
The convergence of frequency versus model size is shown in 
figure 8(b). A logarithmic scale was used along the x-axis to 
accommodate the size of the large finite element models. 

Solution from a commercial code is also included because 
the research-level IFM/ Analyzers software cannot solve an 
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Figure 7. — Stress and displacement solutions for clamped plate, (a) Finite element model for clamped 
plate, (b) Convergence of displacement, (c) Convergence of moment. 




Figure 8. — Frequency analysis of turboprop blade, (a) Finite 
element model, (b) Convergence of frequency. 


TABLE III.— FUNDAMENTAL FREQUENCIES 
BY DIFFERENT METHODS 


Model: 

(number of elements, 
degrees of freedom) 

Frequency, 

Hz 

IFM/IFMD 

Stiffness 

method 

Commercial 

code 

1: (132, 924) 

42.5 

87.5 

87.5 

2: (240, 1620) 

41.3 

75.8 

75.8 

3: (528, 3432) 



55.8 

4: (720, 4752) 



49.5 


eigenvalue problem with 3452 degrees of freedom in the 
model. The brick element of the commercial code is compa- 
rable to that of the stiffness method of the IFM/ Analyzers 
software. The IFM/IFMD code produced a satisfactory 
frequency result for the first model with 924 degrees of 
freedom. (The next model with 1620 degrees of freedom 
yielded only a marginal 3 -percent improvement.) The 
stiffness method and the commercial code produced much 
higher values for the frequency. The commercial code 
produced a frequency value of 87.5 Hz for the first model. 
For the last model with 6864 degrees of freedom the fre- 
quency was 43.1 Hz. Frequency calculated by the stiffness 
method converged with a 4.4 percent residual error as shown 
in figure 8(b). 

The animation of the stress mode shape in IFM/IFMD 
(ref. 22) detected a stress concentration region in the same 
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blade due to a small hole, which is qualitatively shown in 
figure 8(a). Typically, stress mode shape is not available in 
the stiffness formulation. 

3.3 Summary on Structures 

For the three examples given in this paper and from solu- 
tions of many other problems it was observed that the 
IFM/IFMD produced high-fidelity solutions with a modest 
finite element model. The stiffness method required a very 
large model; even then the solution accuracy could not be 
guaranteed. The monopolistic dominance of the stiffness 
method can only be justified because we can solve a very 
large number of equations with a computer. The IFM/IFMD 
methods produced accurate stress, displacement, and fre- 
quency results even for modest finite element models 
because of their simultaneous emphasis on the equilibrium 
equations and the compatibility conditions. The stiffness 
method lacked precision because this method neglected the 
compatibility conditions, especially along the numerous 
boundaries of a finite element model. Even though the 
stiffness method may be preferred at present because of the 
availability of well-written software, the development of a 
commercial code based on IFM/IFMD should be initiated. 

4.0 Design of Structures 

Traditionally the concept developed to design a determi- 
nate structure is also used to design an indeterminate struc- 
ture. Consider for example a determinate truss. The concept 
of fully stressed design (FSD) works very well for a deter- 
minate truss. Typically, the FSD concept is extended to 
design an indeterminate truss. This becomes problematic 
(ref. 13) because an indeterminate truss cannot be fully 
stressed. Rules to design a determinate structure should be 
modified prior to their application to an indeterminate 
structure. The compatibility conditions are key ingredients 
that should be accounted for the modification because the 
CCs imposed implicit constraints in the design of an inde- 
terminate structure. It turns out the CC has considerable 
influence on the fully stressed and fully utilized designs as 
well as in the optimization methods. The implicit constraints 
are intrinsic to the design of an indeterminate structure and 
are independent of the choice of an analysis method used in 
the design calculations. The nature of the constraints is easily 
studied via the IFM, however. The rules thus generated 
should be adopted even when the stiffness method analysis 
tool is used in design calculations. 

Cilley (ref. 26), in 1900 discussed some aspects of design 
of an indeterminate truss. Reinschmidt et al. (ref. 27) 
attempted an analytical explanation through an illustrative 
example. Gallagher and Zienkiewicz (ref. 28) has also 
touched upon the issue. Some of our observations on the 
topic are reported in references 13, 29, and 30. The central 
design issue is illustrated by considering truss examples. 


Extension to indeterminate beams, framework, and other 
types of structures is straightforward, though it would 
require extensive algebraic manipulation. Both traditional 
design and optimization concepts are examined for an 
indeterminate truss, and the discussion is presented under 
three topics: 

(1) Design for strength 

(2) Design for stiffness 

(3) Design optimization 

4.1 Design for Strength 

Design for strength attempts a full stress state for all 
members of the truss. Consider a truss with n bars. For all 
bars in the truss a full stress condition is defined as 

<Ji = g 0 (/ = 1,2, ...,w) (55) 

where a* is stress in the i th bar and a 0 is the strength. A stress 
ratio rule yields the area A t of the z th truss bar as a ratio of 
member force F t to strength: A t = F/g 0 . The rule works well 
for a determinate truss but becomes problematic for an 
indeterminate truss. 

4.1.1 Design of a determinate truss . — Consider a 
determinate truss with n bars and the same number of bar 
areas. The bar areas can be calculated from the governing 
equation of the IFM [S] {F} = {P } by expressing it in areas 
A t = F/g 0 to obtain = {F*}. For a determinate 

structure the governing matrix is equal to the equilibrium 
matrix ([5] = [F]), which is independent of bar areas. In 
other words, a single inversion without iteration yields the 
areas 

{A} = abs {[$*]- V }} (56) 

The formula in equation (56) can be extended to include 
buckling limitations. 

Consider the two-bar determinate truss shown in 
figure 9(a). It is made of steel and is subjected to two load 
components, P x = 10 kip and P y = -20 kip. The bar force is 
linked to the cross-sectional area A through the yield stress 
as F\ = <3 0 A i and F 2 = a 0 A 2 - The IFM governing equation 
degenerates into the equilibrium equations. The EEs are 
modified to obtain the design. 

[*]R1 = {r} => [*]M = =* [B] M =W{P} 

Ivoj 

or {^} = abs|W[ J S]“ 1 {p}| 

H c oJ J 

( 57 ) 
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Figure 9. — Design of indeterminate truss, (a) Determinate, 
(b) Indeterminate. Dimensions are in inches. 

The bar areas are determined from equation (57) in a sin- 
gle step because the equilibrium matrix [ B ] is independent of 
the design variables. The numerical values for the areas are 
as follows: 

j^j tw ° bar _ jo ^ ji n } and weight = 43.4 lbf 

4.1.2 Design of an indeterminate truss . — The truss is 
made indeterminate by adding a third bar as shown in 
figure 9(b). The following design equation is obtained when 
the determinate design concept is extended to the indetermi- 
nate truss: 


[S]{F}={P*} 


or 


{v4} = abs 



(58) 


Equation (58) has to be solved iteratively since the [S] 
matrix is a function of bar areas because of the flexibility 
matrix [G] in the compatibility condition [C][G]{F} = {0}. 


Iterative solution of equation (58) yielded the following 
design: 


( ,) three bar 
{A} =< 


0.71 

0.50 >in. 2 and weight = 43.4 lbf 
0.00 


An extension of the determinate design concept to an 
indeterminate structure returned with a zero area for the third 
bar. Furthermore, design and weight are identical for both 
determinate and indeterminate structures. This observation is 
not coincidental but is common to all indeterminate trusses. 
The implicit compatibility constraint distinguished between 
the designs of the determinate and indeterminate structures. 
This implicit compatibility constraint in deformation, force, 
and stress variables (P) , {F}, and {a}, respectively, is as 
follows: 


Pi - V^p2 +P3 - 0 


AE) X \ AE 


2 + U e ) 3 ° 


(59a) 

(59b) 


a l ~~ a 2 + a 3 = 0 


(59c) 


The compatibility constraint is violated (ref. 13) when a 
fully stressed state (gi = a 2 = a 3 = a 0 ) is attempted: o x - g 2 + 
g 3 => 20 - 20 + 20 * 0. Because of the violation, a bar area 
of zero is obtained for the third bar (A 3 = 0) of the determi- 
nate truss. This observation can be extended for a general 
indeterminate truss with n bars and r redundant members. 
For such a truss ( n , m) there are m = n - r equilibrium 
equations and r = n - m compatibility conditions. A fully 
stressed condition cannot be achieved for an indeterminate 
truss (ft, m). An attempt to fully stress this truss would 
degenerate it to a determinate truss with areas of zero for r 
number of bars: Aj — 0 where j = 1 ,. . ., r = n - m. 

4.1.3 Linking of design variables. — The situation associ- 
ated with bar areas of zero in the design of an indeterminate 
truss can be alleviated by linking the design variables. The 
number of independent variables must be reduced to m for an 
ft-bar truss with r redundant members. For the three-bar truss 
(see fig. 9(b)), we discuss two of several choices to link 
areas: link the areas of bar 1 and bar 3 (A h A 2 , A 3 = A x ); 
likewise link A 2 and A 3 . The designs obtained for the two 
cases are as follows: 


{4 


A\ and A 3 


0.64 

0.66 

0.64 


m. 


and weight = 7 1 .5 lbf 
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Figure 10. — Twenty-bay truss. 


r 0.50 

{A} Al and ^ 3 = < 0.71 > in. 2 and weight = 63.8 lbf 
0.71 

The weight of the structure is increased to 71.5 or 63.8 lbf 
from the weight of 43.4 lbf for the determinate system, when 
areas are linked. Furthermore, the weight of the structure 
depends on the choice of the link strategy. For the three-bar 
truss, weight was increased by 65 percent from the determi- 
nate truss when areas A x and A 3 were linked. The increase 
was 47 percent for the link condition A 2 =A 3 . 

The link strategy has considerable influence on the weight 
of the structure. However, it is not straightforward to deter- 
mine the best link strategy. The bandwidth of the CC pro- 
vides some analytical guidance to develop a link strategy. 
Consider for example the 20-bay truss shown in figure 10. 
This truss has 101 members with 21 redundancies and 
21 CCs. The bandwidth of the CCs are 6 for 20 of them and 
20 for 1 of them. The bandwidth of six pertains to a bay in 
the truss; that is, the first bay has six bars (1 to 6). The areas 
of these six bars are candidates for linking. A simple strategy 
is to link the diagonal bars (A 3 = A 4 ). Likewise, link bars ( Ag 
= A 9 , A 13 = A 14,..., A 98 = A 9 9 ) for other bays. The remaining 
CC pertains to the 20 bars of the bottom chord. A simple 
strategy is to link all such bars ( A 2 = A 7 = A n = ... A 97 ). The 
compatibility bandwidth provided a guideline, but there is no 
unique strategy. It may be tempting to employ linear pro- 
gramming for the purpose. The stiffness method cannot 
provide a link strategy because of the nonavailability of the 
bandwidth of the CCs. 

The CC explained why an indeterminate structure cannot 
be fully stressed. The bandwidth of the CC can be used to 
develop a linking strategy. 

4.2 Design for Stiffness 

Design for stiffness included displacement limitations in 
addition to the strength constraints. The design obtained for 
strength remained adequate if displacement became a pas- 
sive constraint. The design has to be modified when dis- 
placement constraints become infeasible. The design can be 


prorated uniformly to satisfy the displacement limitations. 
The prorated technique can produce an overdesign condition 
because of an underutilization of the strength. The following 
discussion is an attempt to alleviate the overdesign 
condition. 

4.2.1 Stress-displacement relationship . — Displacement 
and stress are dependent upon one another, regardless of the 
bar area of a truss. The dependency relationship (ref. 29) is 
derived for the three-bar truss shown in figure 5. 



(60a) 


(60b) 

o 2 -£(X 2 ) 

(60c) 

03 =~(A, + X 2 ) 

(60d) 

x \ = ~j{ a \ -<*3) 

(60e) 

*2 =-§( 02 ) 

(60f) 


The displacement and stress relationship can be general- 
ized in the following formula: 

Pa 

G k = L C J kX J ( 61 ) 

j=P\ 

The coefficients c jk are independent of the design vari- 
ables, here bar areas. The actual number of displace- 
ment components Xj in equation (61) depends on the 
column bandwidth of the equilibrium matrix [B] that is a 
small number. The maximum bandwidth is four for a 
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two-dimensional truss. It is six for a space truss. In other 
words, a small number of displacements are related to a 
stress component. For the three-bar truss, displacement 
depended on one or two bar stresses. 

One approach is to modify the stress allowable or strength 
for selected members using equation (61). Then proceed 
with the fully stressed design. Let us assume the design for 
strength obtained for the three-bar truss indicated a 
15 percent violation on the limitation specified on the 
displacement component X 2 . From equation (60f), we can 
satisfy the displacement constraint if the corresponding 
allowable stress is reduced by 15 percent (g 2 o = 0.85g 0 ). 
This in turn would increase the area of second bar by 

17.6 percent ( A 2 ™ = 1 . 176^2 ld ) • The design update 

requires no iteration for a determinate truss. Iteration may be 
necessary for an indeterminate truss. The rule is easily 
generalized for a truss with k violated displacements. 

4.2.2 Design formula for displacement limitation. — A 
general formula to design for displacement limitation can be 
obtained by modifying the displacement-force relationship of 
the IFM. 

{X} = [J][G]{F} (62) 

The rectangular matrix [./], is replaced by [[5'] 1 ] with 
appropriate modification to obtain a displacement formula as 
follows: 

|^-| = [[^] r [G]{F} = L]{4 (63) 

Equation (63) is manipulated to obtain an expression for 
the bar areas for displacement limitations as 

[s] {A} = jfLj , or {A} = [s]- 1 {x} (64) 

Equation (64) is combined with equation (58) to obtain 
formulae for stress and displacement limitations: 


obtain a method called modified fully utilized design 
(MFUD). Discussion of the MFUD technique is given in 
reference 3 1 and is not repeated here. However, results for a 
set of nine problems are presented in table IV. Each problem 
is solved twice: first by a design optimization method and 
then by MFUD. 


TABLE IV.— NORMALIZED WEIGHT OBTAINED 
BY SUMT AND MFUD METHODS 


Problem 

Normalized 

weight a 

Number of active constraints 

SUMT b 

MFUD C 

Stress 

Displacement 

SUMT 

MFUD 

SUMT 

MFUD 

3 -bar truss 

1.0 

1.00 

2 

2 

1 

1 

5 -bar truss 

1.0 

1.00 

- 

- 

1 

1 

Tapered 5 -bar 
truss 

1.0 

1.00 

2 

3 

1 

2 

10-bar truss 

1.0 

1.02 

1 

- 

2 

2 

Tapered 10- 
bar truss 

1.0 

1.00 

5 

3 

2 

2 

25 -bar truss 

1.0 

1.00 

2 

4 

4 

4 

20-bay truss 

1.0 

1.02 

13 

13 

2 

2 

60-bar trussed 
ring 

1.0 

1.00 

11 

18 

1 

1 

Geodesic 

dome 

1.0 

1.01 

46 

46 

1 

1 


a Weight is normalized with respect to SUMT solution. 
b SUMT is sequential unconstrained minimization technique. 
C MFUD is modified fully utilized design. 


Consider for example the design optimization of the 
60-bar trussed ring. It was subjected to three load cases. It 
had 60 stress and 6 displacement constraints for each load 
case. The design problem was solved using MFUD and an 
optimization method using sequential unconstrained minimi- 
zation technique (SUMT). The 60 bar areas were linked to 
obtain 25 design variables. For the MFUD method and 
SUMT algorithm, the optimum weights were 308.1 and 
309 lbf, respectively. The active constraint set for MFUD 
method included 1 8 stresses and 1 displacement, against 1 1 
stresses and 1 displacement for SUMT. The MFUD method 
generated a good design solution with few calculations. The 
other eight examples in table IV follow the pattern of the 
ring problem with minor variations. 


14 


stress 


= abs 


V CT 0 ) 


Lf 1 U} 


{4 disp =abs([s] 1 {W}) 
{A} = max({^} stress ,{/l} disp ) 


(65a) 


(65b) 

(65c) 


The formulae in equation (65) are modified for computa- 
tional efficiency and superior convergence characteristics to 


4.3 Design Optimization 

Because of dependent constraints, an optimization algo- 
rithm can encounter difficulty in reaching an optimal solu- 
tion for an indeterminate structural design problem. For 
example, in an n - bar truss, r out of n bar stresses {a} are 
dependent. Stresses are dependent because of the r compati- 
bility conditions, which can be written as j^cjja} ={0} . 

The rxn matrix |^cj is independent of the design variables. 
Likewise stresses and displacements are dependent. The 
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m = n- r number of displacements {X} are dependent on the 
n bar stresses: {x} = [i? {a} . Again the mxn matrix \j3 J 

is independent of the design variables. 

In an optimization algorithm, the calculation of the search 
direction {d} requires the constraint gradient matrix [Vg]. 
An example of the use of the sensitivity matrix to generate 
the direction follows: 

{</} = [0]{Vf} (66) 

r n ['Hvgftvgf [vgf'W 

where [Q\ = . L J — 

0. 5A /{[//]{Vf}} r {[//]{Vf}} 

[^] = M-[vg][[v g f[v g ]]' 1 [v g f 

f is the objective function 

The constraint dependency makes the matrix 

[Vg] [Vg] singular. An optimization algorithm may yield 

a solution despite the singularity condition because most 
often the [Q\ matrix is approximated; it is seldom calculated 
in closed form. It is reinitialized into an identity matrix when 
corruption is suspected. Two approaches have been devised 
to alleviate the singularity issue in optimization. The inde- 
pendent constraints can be segregated (ref. 30), or the design 
sensitivity can be approximated (ref. 32). 

4.3.1 Segregation of independent constraints. — The seg- 
regation of a set of independent constraints from all the 
specified set of constraints is straightforward though it adds 
extra calculations to optimization, which is already computa- 
tionally intensive. The process is illustrated considering the 
three-bar truss as an example, see figure 9(b). The problem 
has three design variables and five constraints. The con- 
straint gradient or the design sensitivity matrix can be 
arranged as 


;{Vai}{Va 2 }{Va 3 }{VX 1 }{VX 2 }; 


reduce to 1" independent 

| dependent 

(67) 

jv^}{vx 2 } 




A lower and upper triangular factorization can separate the 
two independent constraint gradients Vgi and VX 2 , which can 
be used to calculate the search direction. The factorization 
process has to be repeated before the calculation of each 



Number of iterations 


10x10 3 



0 50 100 150 200 


Number of iterations 

Figure 11. — Design optimization of five-bar truss. 

(a) Tapered one-bay truss, (b) Standard optimization. 

(c) Optimization accounting for singularity. 

search direction because the values of the gradients depend 
on the current values of the design variables. 

The singularity issue is illustrated considering the design 
optimization of the five-bar truss shown in figure 11(a). 
Optimum solution is to be obtained for the minimum weight 
condition for stress and displacement constraints. A quad- 
ratic programming algorithm was used to solve the problem. 
The design solutions are graphed in figures 1 1(b) and (c) for 
two cases: 
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(1) Singularity is disregarded (see fig. 1 1(b)). 

(2) Singularity is alleviated through a singular value 
decomposition (see fig. 11(c)). 

Design case 1 , for which independent constraints were not 
separated, exhibited a wide variation in the weight, and it 
failed to converge. Case 2, which used independent con- 
straints, converged to the optimum solution. 

4.3.2 Approximate sensitivity . — Design sensitivity can be 
approximated for an indeterminate structure. The approxima- 
tion alleviates the singularity condition. It also substantially 
reduces the number of calculations in structural optimiza- 
tion. The sensitivity matrix for the stress and displacement 
constraints for an ft-bar truss with r redundant members is 
obtained by differentiating the IFM equations (ref. 33). 
Consider the closed-form sensitivity matrix for stress. It has 
two distinct factors as follows: 




[0] 


[c] 



(70) 


The elements of the diagonal matrix 


■s* 


g 


are given by 



4 % 


(71) 


The diagonal matrices 



and 



are of dimensions nxn , and their elements are 


[Vct]=[Vct 1 ,Vct 2 ,...,Vct„] 


1- 

r . T 

- determinate 

r 

r 

- indeterminate 

- 

1, 

+ 

e-h 

1 i 

1 

A .. 



( 68 ) 


1 F l- 

— , — V , * , and F„ respectively. 

A Af A? Ei 

The determinate factor in the sensitivity of displacement 
can be specialized for an ^-bar indeterminate truss as 


The first factor with superscript “determinate” has to be 
used for both determinate and indeterminate trusses. The 
second factor with superscript “indeterminate” is applicable 
only to indeterminate structures. 

The rank of the nxn sensitivity matrix [Vg] is reduced to 
m = n - r. However, it has full rank n for the first term with 
superscript “determinate.” It is also a diagonal matrix and 
requires a trivial amount of calculation. The second term, 
superscripted “indeterminate,” is responsible for reducing 
the rank, and its generation is computationally intensive. The 
proposition is to use only the first term in equation (68). 

The approximation process can be extended to the dis- 
placement constraint to obtain 


[vx] 



determinate 


+ [[j][G][D]] indeterminate 


(69) 


The proposition is to retain the first term in the displace- 
ment sensitivity. The sensitivity expressions given in equa- 
tions (68) and (69) should be adjusted for the allowable 
limits prior to their use in design optimization. The defini- 
tions of the symbols in equations (68) and (69) are as 
follows: 


[vx] determinate = p i/] ^ dg jj^ 


-FI 

A 2 E 


[A 1 


(72) 


Calculation of the determinate sensitivity for the dis- 
placement essentially requires a back-substitution step with 
the factored form of the inverse of the [S] matrix. The 
inverse matrix calculated to generate force is available at this 
stage. 

Sensitivities of both the stress and displacement contain 
the member forces {F} and the square of bar areas {A}. 
Even for an active displacement constraint, the design in 
essence is modified through the member force. The geomet- 
ric parameters and material properties of the truss are not 
explicitly contained in the sensitivity expression for stress 
because it is a local variable. Because displacement is a 
global variable, its sensitivity expression explicitly contains 
the configuration parameters, the material properties, and the 
design variables. 

The proposed approximate sensitivity of stress and dis- 
placement constraints alleviates the singularity condition 
in the design optimization of an indeterminate truss because 
the coefficient matrix in the determinate factor (see eqs. (68) 
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TABLE V.— OPTIMUM SOLUTION FOR 20-BAY TRUSS DESIGN PROBLEM 
SHOWN IN FIGURE 10 


Sensitivity 

Weight, 

lbf 

Design variables, 
bar area in in. 2 

Active constraints 

Computation time, 
sec 

A, 

A 2 

A 3 

A 4 

A 5 

Stress 

Displacement 

Simple 

2023.2 

3.44 

7.03 

0.41 

1.38 

1.38 

3 

1 

2.5 

Analytical 

2021.8 

3.44 

6.99 

0.41 

1.39 

1.39 

3 

1 

14.8 


and (69)) has a full rank. The optimum solution is reached 
with fewer calculations because the optimization process 
becomes more robust and the sensitivity is generated with a 
trivial amount of computations. The benefit that accrues 
from the use of approximate sensitivity is discussed in 
reference 32. Here, the performance of approximate sensitiv- 
ity is illustrated considering the 20-bay truss shown in figure 
10 as an example. The truss was subjected to three load 
cases. The first load case consisted of forces in the negative 
y-coordinate direction along the bottom chord nodes: 
-40 kip at the midspan and -1 kip at the other nodes. For the 
second load case, all the top chord nodes were subjected to a 
3 -kip force along the x-coordinate direction. For the third 
load case, all the bottom chord nodes were subjected to a 
-3 -kip force along the negative x-coordinate direction. The 
allowable stress was a 0 = 20 ksi. A displacement limitation 
of 0.5 in. was imposed at the midspan nodes 21 and 22 along 
the x- and y-coordinate directions, respectively. The 101 bar 
areas were grouped to obtain five linked design variables. 
The first two design variables A\ and A 2 represented the bar 
area of the top and bottom chord members, respectively. All 
21 vertical bar areas were grouped to obtain the third vari- 
able A 3 . The last two design variables A 4 and A 5 represented 
the bar area of the leading and lagging diagonal members, 
respectively. The five- variable problem had a total of 312 
stress and displacement constraints. Optimum solutions 
obtained by the sequential quadratic programming (SQP) 
algorithm are given in table V. The convergence of weight 
versus the computation time to solution is graphed in fig- 
ure 12. The convergence patterns for both the closed-form 
and approximate sensitivities graphed in figure 12 display 



Figure 12. — Solution time for sequential quadratic 
programming (SQP) algorithm for 20-bay truss. 


similar undulations. The SQP algorithm converged to the 
same solution for the determinate as well as for the analyti- 
cal sensitivities. However, calculations with approximate 
sensitivities converged 590 percent faster than those with the 
closed-form gradients. The approximate sensitivity outper- 
formed the analytical gradients in the design optimization 
problem for the truss. 

4.4 Summary on Design 

The compatibility condition considerably influences the 
design of an indeterminate structure because design is stress 
driven. Ignoring compatibility amounts to a brute-force 
approach that is not likely to produce a robust design that 
industry would use. The infeasibility of fully stressed design, 
singularity in optimization, and sensitivity calculations in 
structural design have been examined for truss-type struc- 
tures. The same three issues should be examined for flexural 
members following the approach that has been established 
for the truss. The behavior study for flexural members is more 
difficult because stress varies along the member length and 
depth. Understanding the role of compatibility will provide an 
insight to behavior of flexural members, leading to a more 
robust design method. The neglected fully stressed (or util- 
ized) design concept should be revived through use of the 
compatibility conditions because it has the potential to become 
an alternate tool for the design optimization method. It is 
simple and is practiced in industry. We have addressed modi- 
fied hilly utilized design (MFUD) in the stress and displace- 
ment constraints for truss design. From solution of many 
design problems, it is observed that properly formulated 
MFUD can match or quite often exceed the performance of 
optimization methods. MFUD should be extended to flexural 
members beginning with continuous beams and then frames 
with a variable depth member. 

5.0 Concluding Remarks 

The completed Beltrami-Michell formulation (CBMF), 
which is a tensorial approach in elasticity (with stress as its 
primary unknown), is presented. The tensorial method places 
simultaneous emphasis on the equilibrium equation and the 
compatibility condition. The displacement method and mixed 
method as well as the total formulation can be obtained from 
this method as specialized cases. Such methods become 
equivalent yielding identical solutions. Fidelity of the stress 
solution can be guaranteed when it is generated via the CBMF 
or methods derived from the CBMF. 
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A displacement method that is not derived from the 
CBMF but formulated independently from a consideration of 
the equilibrium equation should be examined for a compli- 
ance of the compatibility condition in the field and on the 
boundary. 

The observations for the CBMF apply to the integrated 
force method (IFM) for discrete analysis because the internal 
force unknowns being stress resultants retain the seed of a 
tensor. That is, the IFM becomes a tensor type of method 
that emphasizes both stress equilibrium and strain compati- 
bility to yield high-fidelity response even for a modest finite 
element model. The derived dual integrated force method 
(IFMD) displacement formulation, being equivalent to the 
IFM, retained the strength of the IFM. The IFMD should 
replace the popular stiffness method especially because a code 
based on the existing stiffness method can be changed to an 
IFMD code with a small amount of programming effort. 

The equations of the IFM were used to formulate a fully 
utilized design method for stress and displacement con- 
straints. The compatibility condition provides an insight to 
the behavior of an indeterminate structure, like the infeasibil- 
ity of full stress design. It identified singularity in structural 
optimization, simplified the gradient of behavior constraints, 
alleviated singularity in structural optimization, and reduced 
computation time to solution. 


The tensorial approach could not be developed earlier 
because the compatibility condition was not fully understood 
either in elasticity or in structures. The strain formulation 
cannot be ignored because the compatibility concept makes 
solid mechanics a research discipline that is practiced in 
academia and large research organizations throughout the 
world. Without the compatibility concept, the solid mechan- 
ics discipline would degenerate into a few determinate 
analysis courses in applied mathematics. The use of a fully 
developed strain formulation has allowed a free movement 
between different analysis methods and response variables 
(like force (stress) and displacement). High-speed computers 
are extremely helpful but cannot replace the strain formula- 
tion in elasticity and the compatibility conditions in struc- 
tural mechanics. The vacuum created when the researchers 
knowledgeable in strain formulation retire may be difficult to 
fill. The premier academic and research institutions should 
encourage the study of the compatibility condition for both 
linear and nonlinear problems thereby completing the theory 
of the solid mechanics discipline. 

National Aeronautics and Space Administration 
John H. Glenn Research Center at Lewis Field 
Cleveland, Ohio 44135, May 14, 2007 
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Appendix A 
Symbols and Acronyms 


Symbols 


A 

B 

b 

IB] 

[Q 

D 

[D] 

{d} 

E 

F, 

{F} 

f 

fee 

[G] 

G r , G e 

[Vg] 

[J] 

[K] 

l 

h 

h 

M X , My, M X y 

m 

m 

n 

Fc? My, n z 

p 

{P} _ 

P r ,Pe 



q 

R 


m 

r,e 


strain energy 

complementary strain energy 
body force 
equilibrium matrix 
compatibility matrix 
elastic domain 
governing matrix of IFMD 
search direction 
Young’s modulus 
component of {F} 
internal force vector 
objective function 
strain formulation function 
flexibility matrix 

coefficient functions in Green’s theorem 

sensitivity matrix of constraints 

displacement coefficient matrix 

stiffness matrix 

length of boundary 

length of boundary segment 1 

length of boundary segment 2 
moments 

number of displacement unknowns 
mass matrix 

number of force unknowns 
direction cosines of outward normal 
load 

load vector 

external prescribed loads in r and 0 directions, 

respectively 
traction component 
stiffness method load vector 
transverse load 
reaction 

initial deformation vector 
polar coordinates 

number of redundant force parameters 


[Si 

governing matrix of IFM 

T 

temperature 

t 

plate thickness 

u, V 

in-plane displacement 

V 

body force potential 

W 

potential of work done 

w 

transverse displacement 

X 

component of {X} 

m 

displacement vector 

x,y 

Cartesian coordinates 

a 

coefficient of thermal expansion 

(P) 

deformation vector 

IP 0 } 

initial deformation vector 

y 

shear strain 

% 

strain tensor 

M 

material matrix 

n s 

variational functional of IFM 

G 

stress 

Oo 

yield stress 

[Vg] 

sensitivity matrix for stress 

T 

shear stress 

T ij 

stress tensor 

U 

Poisson’s ratio 

9 

stress function 

CD 

circular frequency 

Acronyms 

BCC 

boundary compatibility condition 

CBMF 

completed Beltrami-Michell formulation 

CC 

compatibility condition 

DDR 

deformation-displacement relation 

EE 

equilibrium equation 

FSD 

fully stressed design 

IFM 

integrated force method 

IFMD 

dual integrated force method 

MFUD 

modified fully utilized design 

SQP 

sequential quadratic programming 

SUMT 

sequential unconstrained minimization 
technique 
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Appendix C 

Solution of Elasticity Problem With BCC 


The use of the boundary compatibility condition (BCC) is 
illustrated through the solution of an elasticity problem. In 
addition, the BCCs are listed for a three-dimensional elastic 
continuum. 


The field compatibility condition: 

d O r ! K-ge) _ Q 

dr r 


(C2) 


C.l Use of BCC in CBMF 


The traction equation along the outer free boundary: 


The use of the BCC in the solution of the completed 
Beltrami-Michell formulation (CBMF) is illustrated consid- 
ering the example of the annular plate shown in figure 13. In 
the text solution for the problem was included for the case of 
a thermal load. For simplicity the solution strategy in this 
appendix is illustrated for a mechanical load (P) applied on 
the outer boundary, at (r = a ), as depicted in figure 13. Other 
parameters of the problem remained unchanged. The field 
equations and boundary conditions given in the text are 
adjusted for mechanical load. The traction condition is 
applied at the outer free boundary, while the BCC is imposed 
on the inner boundary, which is fixed. The BCC should not 
be imposed on the outer free boundary. 

The equilibrium equation in the field: 

-T(a r +c>0) = O (Cl) 



<5 r = P air = a (C3) 

The boundary compatibility condition along the fixed 
inner boundary 


Ge - oa r = 0 at r = b (C4) 

The field equations (eqs. (Cl) and (C2)) and the boundary 
conditions (eqs. (C3) and (C4)) yield the stress solution. 
Displacement boundary condition is not used to calculate the 
stress state. Displacement ( u ) back-calculated from the stress 
(or strain) contains an integration constant that is evaluated 
from the prescribed displacement boundary condition 
imposed along the fixed inner boundary. 

u = 0 at r = b (C5) 

The solution is generated in the following seven steps: 

Step 1 : The sum of stress (a r + Ge) is a constant (c 0 ) from 
the field equilibrium equation (Cl): 


G r + G 0 - C 0 

or, a r = Cq — g 0 


(C6) 


Step 2: The field compatibility condition equation (C2) is 
expressed in terms of radial stress o r by eliminating the 
tangential stress g 0 between equations (C2) and (C6) to 
obtain an uncoupled differential equation as 


da r 

dr 


2 

+ — a r 
r 


?0 

r 


(C7) 


Step 3 : Equation (C7) is solved to obtain o r as 


_ c 0 C 1 
6 ~T ~ 

Z r 


(C 8) 


Step 4: Stress o 0 is obtained by back substitution in 
equation (C6) as 


a r 


c 0 C 1 
2 r 2 


(C9) 
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Step 5: The two constants (c 0 and c i) are obtained from 
the traction and boundary compatibility conditions given by 
equations (C3) and (C4), respectively. 


2 a 2 


So_ 

2 


<4 


r r r ^ 
c 0 | C 1 

2 


V 


(CIO) 


= 0 


J 


Solution of equation (CIO) yielded the two constants: 


c 0 = 


2(l + u)< 


c, =■ 


(l + v)a 2 + (l - u)b 2 
(o- \)a 2 b 2 
(l + v)a 2 + (l - v)b 2 


(Cll) 


P 


though one required an integration and determination of the 
constant from the displacement boundary condition, while 
the other bypassed such procedure. The displacement func- 
tion is 


(>-» 2 )p +i> 2 ) ^ 


(l + o)<2 2 +(l-\j)Z? 2 


V r J 


(Cl 4) 


C.2 Boundary Compatibility Conditions in Three 
Dimensions 

For a three-dimensional elastic continuum, the BCCs have 
been derived from the stationary condition of the variational 
functional of the integrated force method (ref. 34). The 
variational calculation yielded three field equations (ref. 35) 
and three boundary compatibility conditions. The BCCs are 
listed in stress, in strain, and in displacement. 


Step 6: The stress solution is obtained from equations 
(C8), (C9), and (Cl 1) as 


[(l + u)r 2 +(l-u)Z? 2 ] f a 2 ^ 


cm = 


(l + u )$ 2 + (l — \S)b 2 
[(l + u)r 2 -(l-u)Z? 2 ] 


Vr 2 y 


(Cl 2) 


(l + yj)a 2 + (l — u)Z? 2 


2 

\r J 


The stress solution is obtained without using the dis- 
placement boundary condition. 

Step 7: Calculation of displacement involves the trans- 
formation of stress to strain. Strain is integrated to obtain 
displacement. The constant of integration is determined from 
the displacement boundary condition. Calculation of dis- 
placement become trivial for the axis-symmetric annular 
plate because of the simple strain-displacement relationship. 


Stains c r and 8e are calculated from stress as 

-u 2 )* 2 ]| 


[(l — ubrMl- 


(l + u )a 2 + (l - v)b 2 


£0 = 


[(l-u 2 ) 


r 2 + (l ■ 


■ o 


)b 2 _ 


(l + u )a 2 + (l - yj)b 2 


2 "\ 
a 

( p ) 

2 

vr y 

UJ 

f 2\ 
a 

' p ' 

2 

vr y 



(Cl 3) 


The displacement can be calculated from the radial strain 

dt/ ( u \ 

z r = — or the tangential strain 8 q =— . Both tech- 

dry V r y 

niques yield the same displacement function, even 


BCC in stress: 


0 

Oz 


\n z (ctj, - ua z - ua x ) - n y (l + u) t . 

+ f~{ n y ( c z “ - uc y ) - n z (1 + *>)' V } = 0 




ua x - uct v 


+— \ nz { a x~ va y~ va ^ 


)-n z ( 1 + o)t 
)-« x ( 1 + u)t zx } = 0 


\ X 5 y ~ UC Z 


0 

¥ 

8 

+ ^-\>lx[Vy-W z -UO x 

BCC in strain: 


{ n y{^x~ 

K( 


(Cl 5a) 


(Cl 5b) 


)-«x( 1 + u Ky} 

)-«y(l + u )v} = 0 


(Cl 5c) 


of 1 \ of 1 

(n z £y--n y Yy Z J +— [n y e z --n zlyz | = 0 (C16a) 


8 ( 1 

n x £ z -~ n zYxz | = ° ( C16b ) 

V 2 


, n z £ r w rYrz + 

8z\ x 2 x x J 8x 


—yi y Z x -- nx y xy J + 


= 0 (Cl 6c) 
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(Cl 7c) 


BCC in displacement: 


d 2 v d 2 w d 2 

n v — v + n z — o ^r^r - 
y dz 2 dy 2 dydz 


(n z v + n y w^ = 0 


d 2 w d 2 u d 2 


n z—T + n x 


dx 2 X dz 2 dzdx 


{n x w+n z u) = 0 


(Cl 7a) 


(Cl 7b) 


0 2 v 

”w +% 


d 2 u 

dy 2 


a 2 

dxdy 


< 


n y v + n x ii 



where o x , a y , a z , T xy , t xz , and z yz are the stress components; s x , 
s y , £z, Yxy, Yxz, and y yz are the strain components; u, v, and w 
are the displacements; and n x , n y , and n z are the direction 
cosines of the outward normal. The BCCs expressed in terms 
of continuous displacement functions are not satisfied 
automatically. 
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